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Abstract. We prove that for a combined system of classical and quantum particles, it is 
possible to write a dynamics for the classical particles that incorporates in a natural way the 
Boltzmann equilibrium population for the quantum subsystem. In addition, these molecular 
dynamics do not need to assume that the electrons immediately follow the nuclear motion (in 
contrast to any adiabatic approach), and do not present problems in the presence of crossing 
points between different potential energy surfaces (conical intersections or spin-crossings). A 
practical application of this molecular dynamics to study the effect of temperature in molecular 
systems presenting (nearly) degenerate states - such as the avoided crossing in the ring-closure 
process of ozone - is presented. 
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1. Introduction 

The possibility of dividing a system of particles into a quantum and a classical subsystem is 
of wide interest in several areas of physics and chemistry. This division is typically, but not 
exclusively, done by considering the electrons to be quantum, and the nuclei to be classical 
(this is the choice that we make in this work, although the observation of quantum effects for 
protons is nowadays a matter of considerable debate [1]). When the gap is large (i.e. when 
the lowest electronic excitation energy is much greater than the highest available nuclear 
vibrational frequency, which depends on the temperature), the division is properly handled 
by making use of the standard Born-Oppenheimer (BO) separation. Ab-initio Molecular 
Dynamics (AIMD) (2) can then be performed for the system of classical nuclei on the ground 
state BO surface (gsBOMD) - and this type of dynamics has been applied countless times in 
the last couple of decades, either directly or by making use of equivalent but more efficient 
schemes such as the Car-Parrinello (CP) technique [3 J or other alternatives fl4j[5]|. These 
dynamics can be used to compute equilibrium averages by assuming ergodicity and computing 
time averages over a number of trajectories. Notice that when the molecular dynamics is used 
to sample the equilibrium distribution, all the dynamics are physically equivalent insofar as 
they produce the correct distribution, and only efficienty criteria may favour one over another. 

If the temperature is of the order of the electronic gap or larger we cannot assume any 
longer that the quantum particles are constantly in the ground state. Indeed, in many physical, 
chemical or biological processes the dynamical effects arising from the presence of low lying 
electronic excited states have to be taken into account. For instance, in situations where the 
Hydrogen bond is weak, different states come close to each other and non-adiabatic proton 
transfer transitions become rather likely at normal temperature J6J. In these circumstances, 
the computation of ensemble averages cannot be based on a model that assumes the nuclei 
moving on the ground state BO surface. 

If transitions to higher energy levels must be allowed, a different type of dynamics must 
be performed. In the realm of first principles MD calculations, two approaches come to 
mind: Ehrenfest dynamics, and surface hopping 0. Their suitability for the calculation 
of equilibrium properties is however still a subject of study Il8l [m4r3l . Also, in density- 
functional theory (DFT), one could perform MD at T ^ 0, working with partial occupation 
numbers to account for the electronic excitations lfl4l . ideally making use of a temperature- 
dependent exchange and correlation functional lfT5Tl . This scheme is however tied to DFT, and 
is hindered by the difficulty of realistically approximating this functional. 

In this work, we first emphasize that the distribution that is often regarded in the context 
of quantum-classical systems as the "correct equilibrium distribution" (p^ 1 ^ in Eq. (IT2]) 
below), despite being commonly obtained through heuristic arguments, is however only a 
zero-th order approximation (in the quantum-classical mass ratio \JmjM) to the canonical 
equilibrium density matrix associated to a rigorous quantum-classical formulation, as shown 
by Nielsen, Kapral and Ciccotti flUEl. Therefore, there is a priori no reason to ask for 
any mixed quantum-classical theory such as, e.g. Ehrenfest dynamics, to yield exactly 
the Boltzmann equilibrium distribution, p^°\ contrarily to what is often required in the 
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literature ||8l[TTHT3l 

However, although 1S onr y an approximation to the correct quantum-classical 

equilibrium ensemble, it is acceptable when y/m/M is small, and the results are then reliable. 
Hence, pj^°^ will also be the target equilibrium distribution in this work. 

Therefore, in what follows, we will write a system of dynamic equations for the classical 
particles such that the equilibrium distribution in the space of classical variables is in fact 
given by Eq. (fT2l) . This is also a goal of surface hopping methods, although it is not fully 
achieved since these methods do not exactly yield this distribution. We will do this by deriving 
a temperature dependent effective potential for the classical variables, which differs from 
the potential energy surfaces (PES) that emerge from BO equations. It is straightforward, 
however, to write an equation that gives the expression for the effective potential in terms 
of these PES. But note that our approach will be based on the assumption that the full 
system of electrons and nuclei is in thermal equilibrium at a given temperature, and not 
on the assumption that electrons immediately follow the nuclear motion (i.e. the adiabatic 
approximation). 

As an example of the interest in going beyond the PES we mention the issue of quantum 
effects in proton transfer lfT6ll . It is a matter of current debate to what extent protons behave 
"quantum-like" in biomolecular systems (e.g. is there any trace of superposition, tunneling 
or entanglement in their behavior?). Recently, McKenzie and coworkers 0~| have carefully 
examined the issue, and concluded that "tunneling well below the barrier only occurs for 
temperatures less than a temperature Tq which is determined by the curvature of the PES 
at the top of the barrier." In consequence, the correct determination of this curvature is of 
paramount importance. 

As we will see, the curvature predicted by our temperature dependent effective potential 
is smaller than the one corresponding to the ground state PES, in the cases in which the 
quantum excited surfaces approach, at the barrier top, the ground state one. Therefore, Tq 
would be smaller than that corresponding to the ground state PES (see Eq. (8) in HI), and 
hence the conclusion in this reference "that quantum tunneling does not play a significant role 
in hydrogen transfer in enzymes" is reinforced by our results. 



2. Method 



We start our discussion with a system of classical particles, which we divide into two 
subgroups, one of them, of mass m, characterized by the conjugate variables (x,p), and 
the other one, of mass M, characterized by the conjugate variables (X,P). If we had only 
one degree of freedom for each subgroup (i.e. only one particle in one dimension), the 
Hamiltonian would read: 

H total (x,p,X,P) = ^ + ^- + V total (x,X) . (1) 
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The generalization of the following to more particles or degrees of freedom is straightforward. 
In the canonical ensemble, the average of any observable 0(x,p,X,P) is computed as: 

{0{x,p,X,P)) = -! dxdpdXdP 0(x,p,X,P)e-^ Hto ^ x ' p ' x ^ , (2) 

where Z = J dxdpdXdPe,~^ Htotsl ^ x,p,x ,p ^ is the partition function. However, if we think of an 
observable A(X,P) that depends only on the variables of one of the particles (let us call it a 
heavy particle, assuming that M ^> m), these two equations can be exactly rewritten as: 

(A(X,P)) = -J dXdPA{X,P)Q-^ft{x,P$) ? ( 3 ) 

with Z = f dXdPe~$ Hetl ( x,p '& . Here, we have defined a temperature-dependent effective 
Hamiltonian: 

H eS (X,P; p) := ~ log J dxdpe^ H ^ x ^ x ^ . (4) 

Therefore, the equilibrium properties of the subsystem formed by the heavy particle can be 
described in a closed manner, with an effective Hamiltonian in which the coordinates of the 
light particle have been integrated out. 

If we now consider the two particles to be quantum, the system will be characterized by 
the canonical operators (x,p,X,P), related by the commutation relations, [X,P] =ih 1 [x jP ] — 
ih, and by a total Hamiltonian H[ 0tSL \{x,p,X,P) whose expression is analogous to Eq. ©, 
except that we now have operators instead of classical variables. 

The key object, as far as equilibrium properties are concerned, is the equilibrium density 
matrix, defined as: 

(3 eq = I e -P" 5 z = Tr[e -p£(i,W)] 5 (5) 

which allows to compute equilibrium averages as: 

(6(x,p,X,P)} = Tr[6(x,p,X,P)p eq ] . (6) 

Like in the classical case discussed before, for observables depending only on the 
operators X,P, the effective Hamiltonian can be defined analogously to Eq. © reading: 



H eS (X,P;P):=-UogTv q 



e -pW totaI (x,p,X,P) 



(7) 



Notice that the integrals in Eq. © are now replaced by the trace operation over the quantum 
Hilbert space denoted by Tr q . 

We are interested, however, in an intermediate case: the heavy particle is classical, 
whereas the light particle is quantum mechanical. For this purpose, it is most suitable to follow 
Kapral, Nielsen and Ciccotti flUEGl and start from the partial Wigner representation ifTTl of 
the full quantum mechanical density matrix, p(f): 

p w (X,P,t):=(2nh)- 1 Jdze' p ^ h (X-z/2\p(t)\X + z/2}, (8) 
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which is an operator only in the light particle Hilbert space, depending on two real numbers 
(X,P). The classical limit can then be taken in a very straightforward manner by considering 
the evolution equation for pw, and retaining only linear terms in y/m/M = /j. The result is [|9l : 

= ~T [Hw, total, PwJ + 2 {{Hw, total, Pw} ~ iPw, Hyy, total} ) • (9) 

In this equation, the {-,•} are the usual Poisson brackets. The Hamiltonian H\y to tai (X,P) is 
the partial Wigner transform of the total Hamiltonian for the full quantum system replacing 
the (X,P) operators by real numbers: 

H w ,toUx,p,X,P) = £- + ^- + y total (i,Z) . (10) 

The equilibrium density matrix in the partial Wigner representation at the classical limit 
for the heavy particle, denoted by p^ should be stationary with respect to Eq. ©. If we use 
this property and expand it: 

CO 

p e ;=l^p e ; (n) , (id 

n=0 



it can then be proved IU01 that the zero-th order term is given by: 

peq (0) _ J_ e -pH H ,, total (x,/3,X, J P) ^ ^ 



with Z = Tr q 



. Note that this object corresponds, at fixed classical 
variables (X,P), with the equilibrium density matrix for the electronic states. However, it is 
only an approximation to the true quantum-classical equilibrium density matrix, since it is not 
a stationary solution to the quantum-classical Liouvillian given in Eq. ©. The error made, 
i.e. the rate of change of the distribution as it evolves in time, is given by: 



dpT ] _ p g 

dt M Z 



Wotal -pfl wtnta1 . --P/Wi ^totaA 
dX + dX J 



1 fdH w .total _Rt7„. Mu 

2 V dX 

Mv 



io dX 



(13) 



It can be seen that it grows with P (quadratic dependence at small (3) and it is proportional 
to the velocity of the classical particle P/M. Therefore, the error becomes smaller as the 
temperature grows. 

With these facts in mind, we will consider Eq. (PT2l) as a reasonable approximation and 
use it, for example, to compute averages of observables: 



(6(x,p,X,P))=Tr q J dXdP6(x,p,X,P)p$ {0) . 



Analogously to the preceeding sections [Cf. Eqs. © and ©], if we are interested in 
observables that depend only on the heavy, classical particle, we can define the temperature 
dependent effective Hamiltonian in the form: 

tf eff (X,P;P) := -llogTrqe-P^^^ . (14) 
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It is then easy to verify that the partition function can be written as: 



Z 



dXdPe-V H ^ x > p '-® , (15) 



and the classical observables can be computed as: 

(A(X,P)) = - J 6XdPA(X,P)e~ ?H ^ p '^ . (16) 

Hence, the quantum subsystem has been "integrated out", and does not appear explicitly in 
the equations any more (of course, it has not disappeared, being hidden in the definition of 
the effective Hamiltonian). In this way, the more complicated quantum-classical calculations 
have been reduced to a simpler classical dynamics with an appropriate effective Hamiltonian, 
which produces the same equilibrium averages of classical observables [Eq. (fT6l) l as the one 
we would obtain using Eq. (fT2l) . and hence incorporates the quantum back-reaction on the 
evolution of the classical variables 

In the case of a molecular system, the total (partially Wigner transformed) Hamiltonian 
reads: 

H total (R 1 P) = T n (P)+H e (R), (17) 

where R denotes collectively all nuclear coordinates, P all nuclear momenta, T„(P) is the total 
nuclear kinetic energy, and H e (R) is the electronic Hamiltonian, that includes the electronic 
kinetic term and all the interactions. The effective Hamiltonian, defined in Eq. (fl4l) in general, 
is in this case given by: 

H eS (R,P;$) = T n (P) - ~logTr q e-f^) := T n (P) + V eff (R; 0) , (18) 

where the last equality is a definition for the effective potential. 

Now, making use of the adiabatic basis, defined as the set of all eigenvectors of the 
electronic Hamiltonian H e : 4(^)1%^)) = E n (R)\ x ¥ n (R)) , we can rewrite V e a(R; |3) as: 



V eff (R;$)=E (R)-^\og 



1 + JV 



(19) 



, -p£»o(«) 

;i>0 

where E n o(R) = E„{R) —Eq{R). This equation permits to see explictly how the ground state 
energy Eq differs from V e ff, and in consequence how a MD based on V e ff is going to differ 
from a gsBOMD. In particular, notice that V e ff(R;$) < Eq(R). Also, note that to the extent 
that nuclei do not have quantum behavior near conical intersections or spin crossings |fT8l , 
nothing prevents us to use this equation also in these cases. 

The definition of the classical, effective Hamiltonian for the nuclear coordinates in 
Eq. ([TBI allows us now to use any of the well-established techniques available for computing 
canonical equilibrium averages in a classical system [jj, given in this case by the convenient 
expression (fT6l) . For example, we could use (classical) Monte Carlo methods, or, if we want 

if. Of course, since H e ff in Eq. ([T8T l depends on T, any Monte Carlo or dynamical method must be performed at 
the same T that H e g was computed in order to produce consistent results 
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to perform MD simulations, we could propagate the stochastic Langevin dynamics associated 
to the Hamiltonian (fT8l : 

Mjkj(t) = -V/Veff(*(f);P) -Mjjkj(t) +MjE(t) , (20) 

where S is a vector of stochastic fluctuations, obeying (E,-(f)) = and (Ei(/i)E/(f2)) = 
2ykgTdjjd(ti — ti) which relates the dissipation strength y and the temperature T to the 
fluctuations (fluctuation-dissipation theorem). 

Indeed, it is well-known that these dynamics are equivalent to the Fokker-Planck 
equation for the probability density W(R,P) in the classical phase space lfl9ll : 

where {-,•} is the classical Poisson bracket. 

Any solution to Eq. (TfTT) approaches at infinite time a distribution W eq (R,P) such that 
dtW^R.P) = 0. This stationary solution is unique and equal to the Gibbs distribution, 
Weq(i?,P) = Z~ l e -P#eff(*^;P) [19]. Thus, the long time solutions of Eq. d2D, and hence 
those of Eq. (l20l) reproduce the canonical averages in Eq. (TT6l) . This property, which is also 
satisfied by other dynamics like the one proposed by Nose ll20~ll2~Ul if the H e ff in Eq. (fT8~l) is 
used, comes out in a very natural way from the present formalism while it is yet unclear of 
other ab initio MD candidates for going beyond gsBOMD f[8 Hm - TT3ll . 



3. Applications 

The most obvious applications of our temperature dependent effective potential and its 
associated dynamics are the processes of proton transfer and thermal isomerization whenever 
low lying electronic excitated states have to be taken into account. These processes are 
ubiquitous in chemistry and biochemistry j6l|24l|25]|. In particular, ab initio molecular 
dynamics of intramolecular proton transfer around conical intersection and excited states is a 
topic of current interest ll26ll and a great tradition 11271 . 

As an example of our approach, we show in the following the difference between our 
temperature dependent effective potential and the gs BO PES for the avoided crossing between 
the open and closed forms of the ozone molecule. We focus on the ring closure of this system, 
as depicted in Fig. Q] (inset). Using the CASSCF method (complete active space self consistent 
field) j[22|. we have computed the PES corresponding to the ground and first excited singlet 
states (l^i and 2 X A\), and the effective potential along the relevant reaction coordinate, the 
ring closure angle (]) - using only those states as they are the only relevant ones in this case. 
At around = § — 83.4°, there is a barrier between two possible minima in the gs PES of 
this molecule. The 2 l A\ electronic state approaches at this point the ground state PES, and in 
consequence one might expect that the effective potential, and its derivatives at the cusp, will 
differ considerably from the gs PES values as the temperature goes up. This can be seen in 
Fig.ffl 
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Figure 1. (color online) PES corresponding to the 1 l A\ (Eq) and 2 1 Ai(£i) states, and effective 
potential at the indicated temperatures, for the ozone molecule. The reaction coordinate is the 
molecular angle. 



The situation shown in this figure is very general. In fact, one can prove from Eq. |T9lthat 
if the first and second derivative of E\q at the barrier top verifies: 



d 2 E 



10 



then: 



d§ 2 
d 2 V eff 



(<M [l+e 



dE 



10 



(to 



3<t> 2 (to;P) 



< 



d 2 E 



dp 



(to; 



(22) 



(23) 



i.e. the curvature by our V e ff is smaller than the gs PES curvature. Note that, in avoided 
crossings (cbo) is approximately zero and therefore the above condition will be verified in 
the most interesting cases. 

Here, we have computed the effective potential energy surface for a given reaction 
coordinate. Once this surface is our disposal, it is a trivial task to perform nuclear dynamics 
on top of this surface. As we mentioned in the previous section, this nuclear dynamics can 
be performed in two ways. First of all, one could perform dynamics for ensembles with the 
Fokker-Planck equation. In this case, one would just use the effective potential that has been 
pre-computed in the equation, and use any of the well-known partial differential equation 
solvers to propagate Fokker-Planck' s equation. Alternatively, one can propagate the nuclear 
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dynamics individually, Eq. (|20l ), and then one would use any of the standard propagators that 
are also well described in the literature. 

However, in many situations it will not be advantageous to pre-compute the effective 
potential energy surface (due to a larger dimensionality, etc), and instead the propagation 
would be peformed "on the fly", i.e. simultaneously to the computation of only those points 
in the surface that are necessary. 

Note that the computational cost for the evaluation of the temperature dependent effective 
potential is necessarily larger than the cost for the computation of the gsBOPES. A number of 
excited state surfaces must be computed, which can be done, for example, with the CASSCF 
technique that we have employed here, or with time-dependent density-functional theory. If 
we wish to perform 'on the fly' MD simulations, we will also need to compute the gradients 
of those surfaces, which can be done with linear response theory - similarly to how it is done 
for the ground state surface. 

4. Conclusions 

In this work we have introduced a temperature dependent effective potential and the 
associated constant-7 dynamics which produce in a natural way the Boltzmann equilibrium 
population of the quantum subsystem and its corresponding back reaction; something pursued 
in recent years by many researchers ll8l [TTI4T31 . We justify our only assumption, the 
equilibrium distribution prescribed by Eq. dT2]) . using the formulation of Nielsen, Kapral and 
Ciccotti flUEEOl. Our approach is particularly useful in the case of conical intersection or 
spin-crossing |[T8l . and does not assume that the electrons or quantum variables immediately 
follow the nuclear motion, in contrast to any adiabatic approach. The fact that, when using our 
effective potential, the height of a barrier in the PES and its curvature near the top decrease at 
avoided crossings makes our work relevant in the context of the transition-state theory ||2~3~1 . 
In particular, this is important if one wants to adequately discriminate possible quantum-like 
behavior of nuclei from simple classical effects due to the direct influence of temperature on 
the potential between two metastable states, for example in biological systems |[Tl[T6ll. 
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